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I. INTRODUCTION 


Nowadays, the inclusive hadroproduction of pairs of isolated prompt photons (diphotons) 
is a subject of intense experimental and theoretical studies. From the experimental point 
of view, this process forms an irreducible background in the searches of heavy neutral res¬ 
onances in the diphoton decay channel, such as Standard Model Higgs boson P and it’s 
Beyond Standard Model counterparts [2]. As for the process itself, it allows us to define the 
set of inclusive differential cross sections over such variables as the invariant mass of the pair 
(M), it’s transverse momentum (pr), azimuthal angle between transverse momenta of the 
photons (A0), rapidity of the photon pair (1' 77 ), Collins-Soper angle in the center of mass 
frame of the photon pair ( 9 ) and a few others [3]. Most of these spectra are measured with 
high precision both at Tevatron [3] and the LHC [4j. 

On the theoretical side, providing the predictions for the above mentioned rich set of 
differential spectra is a challenging task even for the state of the art techniques in perturba¬ 
tive Quantum Chromodynamics (pQCD). While, for the inclusive isolated prompt photon 
production, the p^-spectra from CDF [3], ATLAS [6] and CMS [7] are described within 
experimental uncertainties in the Next to Leading Order (NLO) of conventional Collinear 
Parton Model (CPM) of the QCD |8]. Also, the notably good results where obtained for 
these spectra already in the Leading Order (LO) of Ur-factorization in the Ref. J9J [TO ]. 
In contrast, existing NLO CPM calculations, implemented in the DIPHOX m Monte-Carlo 
event generator, provide very poor description of p t and A0 distributions measured by AT¬ 
LAS [T]. In the CPM, the full NNLO accuracy is required to provide qualitatively reasonable 
description of all distributions m- 

Part of these difficulties can be traced back to the shortcomings of the CPM approxi¬ 
mation, where the transverse momentum of initial state partons is integrated over in the 
Parton Distribution Functions (PDFs), but neglected in the hard scattering part of the 
process. Such treatment is justified for the fully inclusive single scale observables, such as 
deep inelastic scattering structure functions or p^-spectra of single prompt photons and jets, 
where the corrections breaking the collinear factorization are shown to be suppressed by a 
powers of the hard scale [13] . 

For the multi-scale differential observables, there is no obvious reason why the fixed-order 
calculation in the CPM should be a good approximation. Usually the simple picture of fac- 
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torization of the cross section of the hard process into the convolution of hard-scattering 
coefficient and some PDF-like objects is kept, but kinematical approximations are relaxed. 
In the treatment of Initial State Radiation (ISR) corrections in the Soft Collinear Effective 
Theory (SCET) [Tf] or in the Transverse Momentum Dependent (TMD) factorization for¬ 
malism [131 HE] , the transverse-momentum of the initial-state parton is kept unintegrated 
on the kinematical level, but neglected in the hard-scattering part, which is justified e. g. 
when the px of the exclusive final state is much smaller than it’s invariant mass, so that the 
following hierarchy of the light-cone momentum components for the initial-state parton is 
preserved: q T <C |q T | « q ± = x\fS. 

In the opposite limit, when q T <C |qr| ~ q ± = x^/~S, the /cy-factorization [lTj is valid, 
and transverse momentum of the initial state parton can no longer be neglected in the hard 
scattering amplitudes. To obtain the suitable hard scattering matrix element we will apply 
the hypothesis of parton Reggeization, which will be described below. In what follows we 
will refer to the combination of fc^-factorization with hard scattering matrix elements with 
Reggeized partons in the initial state as the Parton Reggeization Approach (PRA). This 
approach is mostly suitable for the study of the production of the final states with high px 
and small invariant mass in the central rapidity region. At high energies \/S' ^ px, such 
final states are produced by the small-x partons, and the resummation of log(l/x)-enhanced 
terms into the unintegrated PDF (unPDF) can be implemented [18] - Clearly, the regions 
of applicability of the TMD and kx- factorization are overlapping, and they should match 
when x —y 1 . 

Turning back to the photon pair production, we can conclude, that neither TMD nor 
&>r-faetorization covers all available range of experimental data. Most of the cross section 
comes from the region where the diphoton has small px and photons fly nearly back-to-back 
in the transverse plane, so additional QCD radiation is kinematically constrained to be soft 
and collinear, and the approach of SCET factorization will be preferable. On the contrary, 
at high px and small A 0 , the A/p-factorization will do a good job, as we will show below. 

The previous attempts to study the prompt diphoton production in /^-factorization [T9], 
El had their own problems. In the Ref. [T9], only LO 2 —> 2 subprocesses where taken into 
account. While, the PRA was used to obtain the gauge invariant expression for the QQ —>- 77 
matrix element with off-shell initial state Reggeized quarks ( Q ), the matrix element for the 
RR —>- 77 with off-shell Reggeized gluons in the initial state was taken the same as in the 
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CPM. In fact, this contribution was overestimated in Ref. [19] due to the erroneous overall 
factor 4 in the partonic cross section of the subprocess gg —y 77 presented in the Ref. [21], 
which have lead to the accidental agreement with the early Tevatron data [22]. This factor 
was carefully checked against the results presented in the literature [23[ El], as well as by 


our independent calculations of the exact RR —> 77 amplitude, described in the Sec. [TV] of 
the present paper. 

In the Ref. [20], the attempt to take into account the NLO 2 —> 3 subprocesses was made, 
but manifestly non gauge invariant matrix elements where used both for 2 2 and 2 —)• 3 

subprocesses. Also, the unavoidable double counting of additional real radiations between 
NLO q*g* —> 577 subprocess and the unPDF was not subtracted, which have lead to the 
questionable conclusion, that no resummation of the effects of soft radiation is needed in 
the small-pr region to describe the data. 

In view of above mentioned shortcomings of the previous calculations, the present study 
has two main goals. The first one is to calculate the real part of NLO corrections to the 
process under consideration in the PRA, and develop the technique of subtraction of double 
counting between real NLO corrections and unPDF in PRA. The second goal is to calculate 
the matrix element of the quark-box subprocess RR —)■ 77 in PRA, taking into account the 
exact dependence on the transverse momenta of initial-state Reggeized gluons. 

The present paper has the following structure, in Sec. |TT] the relevant basics of the PRA 
formalism are outlined and the amplitude for the LO subprocess QQ —>■ 77 is derived. In 


the Sec. Ill the NLO 2 —)• 3 amplitudes are derived and the procedure for the subtraction of 


double counting between NLO real corrections and unPDF is explained. In the Sec. |IV| the 
computation of the amplitude for the quark-box subprocess RR —)• 77 is reviewed, and in 
the Sec. [V]we compare our numerical results with the most recent CDF [3] and ATLAS [4] 
data. Our conclusions are summarized in the Sec. m 


II. BASIC FORMALISM AND LO CONTRIBUTION 

As collinear factorization is based on the property of factorization of collinear singu¬ 
larities in QCD [25], the kx ~factorization is based on the Balitsky-Fadin-Kuraev-Lipatov 
(BFKL) [26j (see (27] [28] for the review) factorization of QCD amplitudes in the Multi-Regge 
Kinematics (MRK), i. e. in the limit of the high scattering energy and fixed momentum 
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transfers. For example, the amplitude for the subprocess g(gi) Tgfe) —■> gfe) +g(qi) + g(gs) 
in the limit when 

S34 —^13? s 45 ^ — ^25) 


where = (g* + qj) 2 , Uj = (g* — g.,) 2 , has the form of the amplitude with exchange of the 
effective Reggeized particle in the f-channel: 


A = 2 S (fi( 93 )7 r ‘f(9.))-^ ( ( — 


1 

1 13 


SO 


^(* 25 ) 


(u(q 5 )Y 2 u(q2)) , ( 1 ) 


where g t i — qi — g 3 , qt2 = gs — g 2 , c, ri, r 2 are the color indices, is the effective qqR 
vertex, r^ | /t r2 (g t i, g t2 ) is the central gluon production vertex iLRg, u;(£) is the gluon Regge 
trajectory. The Slavnov-Taylor identity (g 4 i — gt 2 ) /t Tj;’ l / j. 2 (g t i, g f2 ) = 0 holds for the effective 
production vertex, which ensures the gauge invariance of the amplitude. 

The analogous form for the MRK asymptotics of the amplitude with quark exchange in 
the f-channel was shown to hold in the Leading Logarithmic Approximation (LLA) in [29]. 
For the review of modern status of the quark Reggeization in QCD see Ref. [30]. 

The Regge factor s resums the loop corrections enhanced by the log(s) to all orders in 
strong coupling constant a s , and the dependence on the arbitrary scale s 0 should be canceled 
by the analogous dependence of the effective vertices, taken in all orders of perturbation 
theory. The Reggeized gluon in the t-channel is a scalar particle in the adjoint representation 
of the SU(N C ). In the MRK limit, when all three particles in the final state are highly 
separated in rapidity, the light-cone momentum components carried by the Reggeons in 
f-channels obey the hierarchy qf <C |q t j_| ~ qf, so in the strict MRK limit, the “small” 
light-cone component is usually neglected. 

To go beyond the LLA in log(s), one needs to consider the processes with a few clusters 
of particles in the final state, which are highly separated in rapidity, but keeping the exact 
kinematics within clusters. This is so called Quasi-Multi-Regge Kinematics(QMRK), and 
to obtain the amplitudes in this limit, the gauge invariant effective action for high energy 
processes in QCD was introduced in the Ref. [HI] . Apart from the usual quark and gluon 
fields of QCD, which are supposed to live within a fixed rapidity interval, the fields of 
Reggeized gluons ra and Reggeized quarks [32] are introduced to communicate between 
the different rapidity intervals. To keep the t-channel factorized form of the amplitudes in 
the QMRK limit, the Reggeon fields have to be gauge invariant, which leads to the specific 
form of their nonlocal interaction with the usual QCD fields, containing the Wilson lines. 
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Gauge invariance of Reggeon fields also ensures the gauge invariance of the effective emission 
vertices, which describe the production of particles within the given interval of rapidity, as 
it was the case in (JT|). The Feynman Rules (FRs) of the effective theory are collected in the 
refs. [[32] 133], but for the reader’s convenience, we also list the FRs relevant for the purposes 
of the present study in the Fig. [l] To compute the hard-scattering matrix elements in PRA 
one have to combine the FRs of the Fig. [l]with the usual FRs of QCD and QED, and use 
the factors for the Reggeons in the initial-state of the hard subprocess, also defined in the 
Fig.0 to be compatible with the normalization of the unPDF described below. 

Recently, the new scheme to obtain gauge-invariant matrix elements for k t factorization, 
by exploiting the spinor-helicity representation and recursion relations for the tree-level 
amplitudes was introduced [3T, 35|. This technique is equivalent to the PRA for the tree- 
level amplitudes without internal Reggeon propagators, however, the construction of the 


subtraction terms in the Sec. Ill requires the usage of the FRs of the refs. 

So far as, the form of the central production vertex, propagators of Reggeized gluons 
and Regge trajectories, depends only on the quantum numbers of the Reggeon in the t- 
channel, and do not depend on what particles are in the initial state, the cross section of 
the production of particles to the central rapidity region in the inelastic pp collisions can be 
written in the form 53: 


h3 


Xi .1 7T J X 2 J VT 


$j(x 2 ,t 2 , p 2 F )dcr ij (q 1 ,q 2 ), ( 2 ) 


where the sums are taken over the parton species, qy 2 = X\$Pi ,2 + are the momenta 
of the partons, q 2 2 = — qf - 12 = — P. 2 , -Pi ,2 are the four-momenta of the protons, 2P\P 2 = S, 
da is the partonic cross section with Reggeized partons in the initial state. In what follows 
we will often use the Sudakov decomposition of the momenta: 


k — - (n + k + n fc + ) + kr, 

where k± = k° ± k 3 , Pi = \fSn~/ 2 , P 2 = \fSn + / 2 , (n + ) 2 = ( n ~) 2 = 0, n + n~ = 2. 

The unPDF <hj(x, q^, p 2 F ) is unintegrated over the transverse momentum q|,, but still 
integrated over the ’’small” light-cone component of momentum, so this light-cone compo¬ 
nent is neglected in the hard scattering part, which is hence formally in the QMRK with 
the ISR and therefore it is gauge invariant. The exact kinematics will be restored by the 
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higher-order QMRK corrections. The factorization scale /j ,2 F is introduced to keep track of 
the position of the hard process on the axis of rapidity. 

The unPDF is normalized on the CPM number density PDF via: 

J dt$i(x, t, h 2 f ) = xfi(x, n 2 F ). 

In the case of an inelastic scattering of objects with intrinsic hard scale, such as photons 
with high center-of-mass energy and virtuality, the evolution of the unPDFs is governed by 
the large log(l/x) and they satisfy the BFKL evolution equation [26]. In proton-proton colli¬ 
sions, the initial state does not provide us with intrinsic hard scale, therefore, the fey-ordered 
Dokshitzer-Gribov-Lipatov-Altarclli-Parisi (DGLAP) [36] evolution at small fey should be 
merged with rapidity-ordered BFKL evolution at high-fey final steps of the ISR cascade. 

The last problem is highly nontrivial and equivalent to the complete resummation of the 
log(fey )-enchanced terms in the BFKL kernel. A few phenomenological schemes to com¬ 
pute unPDFs of a proton where proposed, such as the Ciafaloni-Catani-Fiorani-Marchesini 
(CCFM) approach [18], the Bliimlein approach [37] and the Kimber-Martin-Ryskin ap¬ 
proach [38J. In the LO calculations in PRA, the definition of the hard-scattering coefficient 
da is independent on the approximations made in the unPDF, so any unPDF can be used, 
and spread between them gives the theoretical uncertainty. In fact, the recent studies [39, TD] 
show, that in the realistic kinematical conditions, the LO calculations with KMR and re¬ 
cent version of the CCFM unPDFs ra give very close results. At the NLO, one should 
develop the proper matching scheme between the unPDF and corrections included into the 
hard-scattering kernel, which introduces a difference in treatment of different unPDFs. 

In the present paper we will work with the version of the KMR formula for the unPDFs, 
described in the Ref. [12] • The KMR prescription introduces the simplest possible scenario, 
where the fey-ordered DGLAP chain of the emissions is followed by exactly one emission, 
ordered in rapidity with the particles produced in the hard subprocess. Due to the strong 
fey-ordering of the DGLAP evolution, the transverse momentum of the parton in the initial- 
state of the hard subprocess is approximated to come completely from the last step of the 
evolution. With this approximations, one can obtain the unPDF from the conventional 
collinear PDF as follows: 
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<f>i(x,q^,n 2 ) = ^2 J dzTi(q 2 , p 2 ) ~ Y^Pij(z)fj (yW 2 ) 9 ( A * Mt ,^ 2 ) - -) , (3) 

rr ^ 

where Pij{z)- DGLAP splitting function, g 2 = g|,/(l — z) - virtuality of the parton in the 
t-channel, the Sudakov formfactor T t is defined as: 


Ti(q 2 , p 2 ) = exp { - J ^^4^ Y1 J d ^ p ij(0 9 { A ij(k 2 {l - £),p 2 ) ~ f) 1 , (4) 

q 2 l ’ j o J 

and the ordering in rapidity between the last parton emission and the particles produced in 
the hard subprocess is implemented via the following infrared cutoff 




h 


d jj "P (1 djj ) 


P + qr 

In the present study, we use the version of the KMR formula ([3]) with LO DGLAP splitting 
functions, but NLO PDFs as a collinear input, because, as it was shown in the Ref. m the 
usage of the NLO PDFs and the exact scale q 2 are the most numerically important effects 
distinguishing the LO KMR distribution of the Ref. [38j and the NLO prescription of the 
Ref. [02] • Also, as it will be shown in the Sec. Ill the usage of the LO DGLAP splitting 


functions is compatible with the PRA, while at the NLO, the splitting functions should be 
recalculated using the effective theory of the refs. [211 [32]. 

The effects of the Sudakov resummation are known to be dominant in the doubly- 
asymptotic region, q? p and z -)■ A, which is most important in pp collisions, therefore 
we use the Sudakov form-factor in Q. The opposite limit, when z —» 0 and qx <C p, is 
captured by the Regge factor but the proper procedure of matching of the double- 

logarithmic corrections, between Sudakov and Regge factors is also beyond the scope of the 
present study. 

Now we are at step to discuss the LO and NLO contributions to the prompt photon pair 
production in the PRA. There is only one LO ( 0(a 2 a °)) subprocess: 


Q ( q i) + Q(q 2) ->• 7(93) + 7(94) 


(5) 


The set of Feynman diagrams for this subprocess is presented on the Fig. [2j The amplitude 
of the process ([5]) obeys the Ward identity of Quantum Electrodynamics (QED), and the 






amplitude squared and averaged over the spin and color quantum numbers of the initial 
state, which was obtained in a first time in the Ref. [12], has the form: 


\A(QQ —► 77 ) | 2 = ^ 7 r 2 e 4 a 2 - (w 0 + w A S + w 2 S 2 + w 3 S 3 }, (6) 

where a 3 = q 3 /VS, a 4 = qt/y/S, b 3 = q 3 /VS, b A = q A /VS, s = (91 + g 2 ) 2 , t = (q 1 - q 3 ) 2 , 
u = (qi — q 4 ) 2 , x 4 = a 3 + a 4 , x 2 = & 3 + b A , cc = e 2 /(47 t), e ? is the electric charge of the quark 
in the units of the electron charge and the coefficients Wi can be represented as follows: 

wo = tit 2 (ti + t 2 ) - tu(t + u), 

— Wi — tit 2 (a 3 - a 4 )(b 3 - b A ) + t 2 x A (b A t + b 3 u) + 

T t\X 2 (^a 3 t T a A u) T tu(a 3 b 3 T 2 a A b 3 T c 2a 3 b A T o 4 & 4 ), 


- w 2 = b 3 b A x\t 2 + a 3 a A x\ti + a 3 b A t(x 2 a 3 + a A b A ) + a A b 3 u(a 3 b 3 + a A x 2 ), 

-w 3 = a 3 a A b 3 b 4 (a 3 b A + a A b 3 (^j ). 

Taking the collinear limit t lj2 —^ 0 of (|6]) , one can reproduce the standard CPM result for 
the amplitude 99 —> 77: 

\AQQ 77)I 2 = y7r 2 eja 2 ^ . 

In the next section we will discuss the tree-level NLO corrections. 

III. REAL NLO CORRECTIONS 

The tree-level NLO ( 0 (a 2 a ])) subprocesses are: 

Q{qi) + R{q2) -a 7(93) + 7(94) + 9(95), (7) 

Q{Ql) + 0(92) -:> 7(93) + 7(94) + 9(95)- ( 8 ) 

The sets of Feynman diagrams for them are presented in the Figs. [ 3 ] and [ 4 ] The FRs 
of the Fig. [l] where implemented as the model hie for the FeynArts m, Mathematica 
based package, and the computation of the squared matrix elements was performed using 
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FeynArts, FeynCalc [36], and FORM programs. It was checked analytically, that the ampli¬ 
tudes for the NLO subprocesses ([7]) and Q obey the Ward (Slavnov-Taylor) identities with 
respect to all final state photons (gluons) independently on the transverse momentum of the 
initial state Reggeized partons. Unfortunately, the obtained expressions are too large and 
non-informative to present them here. 

The squared matrix element of the subprocess ([T]) contains the collinear singularity, when 
the three-momentum of the quark becomes collinear to the three-momentum of one of the 
photons. This collinear singularity can be absorbed into the nonperturbative parton-to- 
photon fragmentation function, and then, the theoretical cross section is represented as the 
sum of direct contribution, where the collinear singularity is subtracted, according to e. g. 
MS scheme, and fragmentation contribution, which is equal to the convolution of the cross 
section of the parton production in pQCD and the parton-to-photon fragmentation function. 
Experimental (hard-cone) isolation condition require the amount of hadronic energy within 
the photon isolation cone of the radius R to be smaller than the fixed value E^ so ' 1 0 ( 1 ) 
GeV: 


E^ ad) (r < R) < E ( t ISO) , (9) 

where r = \JAg 2 + A0 2 is the distance in the pseudorapidity-azimuthal angle plane, 
E^ ad \r < R) is the amount of the hadronic transverse energy within the isolation cone 
around the photon. This isolation condition strongly suppresses the fragmentation compo¬ 
nent, but at high energies, fragmentation is still non-negligible, constituting up to the 20 % 
of the cross section HU- 

The proper treatment of the collinear singularity, considerably complicates the analytical 
computations both in the NNLO of CPM [12] 07], and in the NLO of PRA. Since in the PRA, 
the part of transverse momentum is provided to the hard subprocess by the phenomenological 
unPDFs. To avoid this difficulties, one can define the direct part of the cross section in the 
infrared-safe way, using the smooth-cone isolation condition [38]: 


r < R 


- ](had ) 


r) < E, 


(ISO) 


\{r), 


( 10 ) 


where x( r ) = ( i- cos( j?!) ) > n — 1/2- The isolation condition (10) is easy to implement into 
the process of Monte-Carlo integration, and it makes the cross section of the subprocess 
(j7j) finite, because the collinear singularities associated with the initial-state are regularized 
by the unPDF. Applying the smooth-cone isolation to our calculation we are completely 
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eliminating the need in the fragmentation component [471 08], but of course this isolation 
do not match to the experimental one. However, as it was shown in the Ref. m, the cross 


section obtained with the isolation condition (10) is a lower estimate for the direct plus 
fragmentation cross section, obtained with the hard-cone isolation. Numerically, for n = 1 
this estimate is very good, since it reproduces the NLO results with standard isolation with 
the accuracy of 0(1%) [07] . Having in mind, that we are going to discuss 0(50% — 100%) 


NLO effects, we will apply the isolation condition (10) to our present calculations. 

The cross section of the subprocess ([8| is also finite, because the Sudakov formfactor 


decreases in the region q^ <C /r 2 faster than any positive power of Qt and therefore regularizes 
the collinear and soft singularities of the matrix element of the subprocess ([8]) in the limit 
of qT5 -t 0. 


In the factorization formula ([2]), the part of the ISR, highly separated in rapidity from 
the particles produced in the hard subprocess is included into unPDFs, and the effects of the 
additional radiations close in rapidity to the hard subprocess, should be taken order by or¬ 
der in a s in the hard-scattering coefficient. Therefore, the corresponding MRK asymptotics 
should be subtracted from the NLO QMRK contributions ([7]), and (J8]) to avoid the double 
counting, when the additional parton is highly separated in rapidity from the central region. 
The analogous procedure of the “localization in rapidity” of the QMRK contributions was 
proposed in the refs. [43, 44]. To be compatible with our definition of the KMR unPDF (|3]) , 
this subtraction term should interpolate smoothly between the strict MRK limit, when addi¬ 
tional parton goes deeply forward or backward in rapidity with fixed transverse-momentum, 
and collinear factorization limit, when the initial-state partons are nearly on-shell and ad¬ 
ditional parton has a small transverse momentum but it’s rapidity is arbitrary. Below, such 
subtraction term is constructed in close analogy with the High Energy Jets approach [49]. 

The Feynman diagrams for the subtraction terms, required for the squared amplitudes 
of the subprocesses ([7]) , and ([8]) are shown in the Fig. [5] and can be easily written according 
to the FRs of the Fig. [l] To extend the applicability of the subtraction terms outside of the 
strict MRK limit, one have to implement the exact 2 4 3 kinematics for the subtraction 
terms, taking into account the exact f-channel momentum in the propagator of the Reggeized 
quark. In what follows, we will refer to the amplitudes with the Reggeon propagators and 
vertices, but without kinematical approximations, characteristic for the MRK, as modified 
MRK (mMRK) amplitudes. As it was checked explicitly, the implementation of the exact 
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kinematics do not destroy the gauge invariance of the subtraction terms with the Reggeized 
quarks in the {-channels, presented in the Fig. [ 5 } as it was the case for the rnMRK amplitudes 
with the Reggeized gluons in the t-channels in the Ref. [49]. 

The last ambiguity, which we have to fix in the definition of our mMRK amplitudes is 
the position of the -projector in the numerator of the propagator of the Reggeized quark. 
In the MRK limit, the “small” light-cone component of the Reggeon momentum can be 
neglected and the projectors P± commute with q t under the sign of the trace, but outside 
of this limit, this is not true anymore. To fix this ambiguity, let’s study the amplitudes of 
the mMRK subprocesses in the Fig. [6] Explicitly, they have the forms: 

■A-fjgQ -> qg) = g 2 s v(ql) (V - T b P + \ ( 7 ^ - qt~^] T a v(q 3 ), 

■A^MQ 99) = 9 2 s v(ql) (ju ~ Q 2 T b P + \ + qj^ T a u(qx), 

\ qt Qa J Qt V Qt J 

— n / 77,+ fi~ \ || 

A^QQ 9 ) = g s v(qb (lu ~ - Qt— ) T a u(q\). 

V Qa Qa j 

Taking the squared modulus of these amplitudes and averaging them over the initial-state 
spin and color quantum numbers, we get: 


(11) 

( 12 ) 
(13) 


I A(gQ -> qg )| 2 = g. 


s +f 2 2P 


qg\ 


I AqQ -»■ gg )\ 2 = 9. 


s + t 2 + t zt 
2 S + t 2 2 P qq (z) 


s + t 2 + t 
CaC f 


zt 


-\A(QQ —t g )| 2 , 


\A(QQ^g)\\ 


\A(QQ -)■ g)\ 2 — g 2 s ^ 2 (q*j_ + ^ 2 ), 


(14) 

(15) 

(16) 


where we have taken the limit q\ = 0 to facilitate the study of the collincar singularity, 
qi = —1 2 , the invariants s, t, u are defined after the Eq. (@h qa. = (qT2-qT 4 ), 2 = 1 ~ot/of 
and P qg (z) = l(z 2 + (1 — z) 2 ), P qq (z) = are the LO DGLAP splitting functions. 

When z<l and f-fixed, the partons 3 and 4 are in the MRK. In the opposite (collincar) 


limit t —> 0, the squared amplitudes (14), and (15) factorize into the collinear singularity 


with the corresponding DGLAP splitting function and the squared amplitude (16). From 
this example one can conclude, that the factor q t should be taken together with the vertex 
of the MRK-emission to correctly reproduce the collinear singularity of the amplitude. This 
prescription is denoted by the crosses on the quark propagators in the Figs. EH and [6] 


The squared amplitudes (14, 15) can also be used to explain the structure of the fac¬ 
torization formula (J2| and the unPDF (|3]). The presence of the exact DGLAP splitting 
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functions in ( p~4| ) , and (15) corresponds to the usage of the exact splitting functions in the 
nnPDF ([ 3 ]). Factor z in the denominators of ( |l4| [TH] ) is nothing but a flux factor of the t- 
channel partons, which tells us, that for the Reggeized partons one should use the same flux 
factor I = 2 Sx\X2 as for the CPM partons. After the integration over the “small” light-cone 
component in the definition of the unPDF, the additional factor 1/(1 — z) appears, which 
converts t into t(l — z) = q)y, that’s why the q f* and not q 2 stands in the denominator of 

0- 


The rapidity ordering conditions are imposed in the subtraction terms for the subprocess 
(| 8 ]) (see Fig. [5j lower panel), while for the case of the subprocess ([7]) the rapidity of the 
quark in the final state is unconstrained (Fig. [5j upper panel). This corresponds to the 
fact, that in the KMR unPDF, the radiation of the gluon is ordered in rapidity with the 
particles, produced in the hard subprocess, while for the quark it is not the case. So, the 
mMRK terms constructed according to the Feynman diagrams in the Fig. [5] are completely 
well-defined and correspond to the definition of the unPDF (|3]). 


IV. THE QUARK-BOX CONTRIBUTION 


The subprocess: 


R(q 1) + r ( q 2) -a 7 (? 3 ) + 7(94), 


(17) 


is described by the quark-box amplitude, and is formally NNLO (0(a 2 a 2 )), but it it’s 
contribution to the total cross section is expected to be comparable with NLO contributions, 
due to the enhancement by two gluon unPDFs. The helicity amplitudes for the subprocess 
0 could be written as: 


A{RR, A 3 A 4 ) = //=«.,«< 2 <,(A 3 )7 4 (-A 4 )V<'“««'“, 


(18) 


where A 3 , A 4 are the hclicities of the photons in the final state and the fourth-rank vacuum 
polarization tensor has the form: 


1M2M3M4 


d 4 q 


(?3 ++ 94 


v. (<?-9 i) 2 (9 + 92 -94) 2 (g + 92)V 

, /I3 -H- /I4) + (q 2 <-> —94, M2 V)- /i4)} , 


(19) 
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where the factor 2 takes into account the diagrams with the opposite direction of the fermion 
number flow. The following overall factor is taken out from the amplitude (18): 


e 2 gl s ab ( v 2 \ 

( 2vr ) 4 2 \ 2 f «) ■ 

We take the polarization vectors for the final-state photons in the form: 


where 


^(A) 


1 

71 


« + iH) 


5 


K = ~r ((<?3<?4)g2 - 3 - fegsK), 


A 


1 

A' 


__ e M2 9394 


and A = \Jstu — stit 2 /2. 

In the Ref. [52] we managed to obtain the compact result for the helicity amplitudes of 
the process 7 R —>■ 7 g, and explicitly demonstrate the cancellation of the spurious collincar 
singularity 1 / (^ 1 ^ 2 ) in the squared amplitude. For the process it turns out to be im¬ 
possible to obtain the reasonably compact results, and the task is actually to obtain the 
representation for the helicity amplitudes which will be feasible for the numerical evaluation 
at all. To do this, we observe, that exploiting the Ward identity gfif Af w/t2A13M4 = 0 for the 
tensor (19), one can make the following substitutions in (pT|: 


Q 1 Q 2 
4 a /ht‘ 


: (n ) M 1 (n + ) M2 —> n^nf 12 


T 2 ’ 


where n T 1,2 = qn, 2 /\/h~ 2 - 

To get rid of the e-tensors and directly pass to the Passarino-Veltman reduction for the 
Feynman integrals with the scalar products in the numerator, we exploit the same trick as 
in Ref. [52] . We decompose the four-vectors 171,2 as follows: 


n T i = /7 )( ?i + / 7 )( ?3 + / 7 )( ?4 + 7i n y, (20) 

«T 2 = Po ) q2 + A3 2)( ?3 + P?qA + 72^2/, (21) 

and the vector n y is introduced via it’s scalar products: n 2 y = —1, n y q 2 = n y qs = n y q 4 = 0 . 
The coefficients of this decomposition can be straightforwardly expressed through the 
Mandelstam invariants, transverse momenta of particles and azimuthal angles. After the 
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Passarino-Veltman reduction, the helicity amplitudes where represented as a linear combi¬ 
nations of two, three and four-point scalar one loop integrals, and the cancellation of the 
Ultra-Violet (UV) and Infra-Red (IR) divergences was checked both analytically and nu¬ 
merically. The coefficients of this decomposition depends on 5 invariants s, t, u, t±, t-z and 
8 coefficients /3^\ 7 ^ 2 , he. 13 parameters in total. They can be represented as rational 
functions with tens of thousands terms in the numerators. It turns out, that just to reliably 
check the cancellation of the UV and IR divergences numerically, one have to compute this 
coefficients with 30 digits of accuracy at least. 

Also, it was checked, both analytically and numerically, that the collinear limit for the 


squared helicity amplitudes (18), defined as: 


2n 

[ ~TY~W *li m J^RR, A3A4)| 2 = - ^ |v 4 .cpm(AiA 2 , A 3 A 4 )| 2 , 

J \IttY t i, 2 ->o 4 z — J 

0 Ai, 2 =± 

holds, where the |« 4 cpm(AiA 2 , A 3 A 4 )| j is the squared helicity amplitude of the process gg —)• 
77 in the CPM. The numerical check of the collinear limit was performed by the technique 
described in the Ref. [ 52 ]. The numerical results for the subprocess will be presented 
in the next section, and the FORTRAN code for the calculation of the helicity amplitudes and 
differential cross sections of the process 0 , as well as for the 2 —> 2 subprocess (| 5 ]) and 
2—^3 subprocesses ([ 7 ]), and (J 8 ]) is available from authors on request. 


V. NUMERICAL RESULTS 


The differential cross section of the 2 —)■ 2 subprocesses can be represented as 

follows: 


da 

dq T2 ,dq TA dA(pdy^dy A 


2?r 

J dt\ j d(j)i E $i(xi,ti,A*F)%(^ 2 ,i 2 ,A*p) 

0 


qT3qr4\Aij\ 2 

2(27r) 3 (SxiX2) 2 ’ 


( 22 ) 


where the factor 1 / 2 ! takes into account the identical nature of the photons, qn = qri|, Vi 
are the rapidities of the hnal-state particles, Acft is the azimuthal angle between transverse 
momenta of the photons, (pi is the azimuthal angle between qp 3 and qn, f 2 = (qp 3 + qp 4 — 
qn) 2 , %i ,2 — + qT A e ±V4 )/\/S. The spectra differential in the diphoton invariant 

mass (M) and diphoton transverse momentum p T = \Jq^ z + q^ 4 , could be obtained using 
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the following substitutions: 


dM = c ^ir dqT4 ’ 

dp T = —dq Ti) 

Pt 


(23) 

(24) 


where C = cosh(y 3 — y 4 ) — cos(A</>), D = \Jp\ — g^ 3 sin 2 (A 0 ), qxA = M 2 /(2qT 3 C) for the 


case of (23) and qx 4 — D — qx 3 cos(A(j)) for the (24). 

For the 2 —)• 3 subprocesses ([7| | 8 b , the formula for the differential cross section reads: 


da 1 f 

dq T3 dq T 4dA(j)dy 3 dy4dy 5 2 !,/ 1 


2t r 


2 t r 


d(j> i / dd ‘2 / dcf) 2 x 


x ]r *,(*„ t„ i 2 , uo 16 p^^' 2) 

l 3 


(25) 


and the differential cross sections over the y-r and M could be obtained using the substitu¬ 


tions (23), and (24) as in the 2 —y 2 case. 

For the numerical computations we use the modified KMR unPDF (J3]) with the MSTW- 
2008 NLO PDFs [50] as the collinear input. Also, the value for the fine-structure constant 
a = 1/137.036 was used in the calculations together with the NLO formula for the a s 
with a s (M z ) = 0.12018 and flavor thresholds at m c = 1.4 GeV and = 4.75 GeV. The 
choice of the factorization and renormalization scales /j,r = fip = C.M commonly used in 
the literature 0 m HB Q2] was adopted, where the default value for £ = 1 and the values 
£ = 2 ±x where used to estimate the scale uncertainty of the calculation, which is indicated 
in the Figs. [ 8 ] — 12 as the gray band. The numerical computations where performed mostly 
using the Suave adaptive Monte-Carlo integration algorithm with the cross-checks against 
the results of Vegas and Divonne algorithms implemented in the CUBA library [5TJ . 

Before the presentation of the comparison of our calculations with experimental data, 
let us discuss the contribution of 2 —> 3 subprocesses (7|8) after the subtraction of double 
counting, discussed in the Sec. |III| In the left panel of the Fig. [7] the y^-spectra for the NLO 
contributions 00 in the ATLAS-2013 kinematical conditions (see the second column of 
the table |T|, is presented together with the corresponding mMRK subtraction contribution. 
For the CDF-2012 kinematics of the Ref. [3] (first column of the table [T| the qualitative 
picture is the same. 

From the Fig. [7] one can observe, that mMRK subtraction term reproduces the exact 
contribution of the NLO subprocess (| 8 |) with the 10% accuracy and constitutes more than 
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pp, CDF-2012 [3] 

pp, ATLAS-2013 g] 

\/S = 1960 GeV 

VS = 7000 GeV 

?T3,4 15, 17 GeV 

<?T3,4 > 22, 25 GeV 

1 2/3,4) < 1-0 

12 / 3 , 41 < 1.37, 1.52 < 1 2/3,4 1 < 2.37 

R = 0.4, E { t SO) = 2 GeV 

R = 0.4, E^ SO) = 4 GeV 


TABLE I: Kinematical conditions for the CDF and ATLAS datasets. 


50% of the cross section of the subprocess ([7]) for the pr > 50 GeV. As the right panel of the 
Fig.0 shows, for the subprocess ([?]) the significant deviation from the mMRK asymptotics 
starts only for Ay = y 5 — V ,, 7 < 2.0 while for the larger values of Ay, the the QMRK 2 —)■ 3 
cross section is well described by the mMRK asymptotics. For the subprocess (J8]) the 2 —)• 3 
cross section is reproduced by the mMRK approximation for all values of Ay. Consequently, 
more than 50% of the cross section of the subprocess ([7]) and almost all contribution of the 
subprocess ([ 8 ]) will be canceled by the subtraction term. Having this in mind we do not 
include the contribution of the subprocess (| 8 ]) in the further calculations. 

The squared amplitude for the subprocess (J7]) can be safely integrated from q T5 = 0 in 
(25). The cross section for the subprocess (( 81 ) is also finite, but for the small values of qTs, 


most of the cross section is accumulated at tip ~ 1 GeV 2 , which is nothing else than the 
manifestation of the usual infrared singularity for the radiation of the soft gluon. For this 
reason the cutoff qr 5 > 2 GeV was imposed to produce the lower panel of the Fig. [7| The 
dependence of the cross section on the small-f behavior of the unPDF is unphysical and will 
be canceled away by the NLO real-virtual interference contribution. 

Now we are in a position to compare the predictions of our model with the experimental 
data of the refs. pm- For the comparison, we choose three main observables, clcr/dpT, 
da/dAcp, as ones, where the fixed-order calculations in the CPM are experiencing the great¬ 
est difficulties and da/dM as the benchmark CPM observable, for which the multi-scale 
nature of the process under consideration is less important. The comparison for the other 
observables will be discussed elsewhere. 

In the Fig. [ 8 ] the p^-spectra of the photon pair, measured by the CDF Collaboration is 
presented. For this dataset, three data samples are provided, the inclusive one and two data 


17 








samples with the additional kinematical constraint px < M or px > M imposed. One can 
note, that the inclusive data and data for the px < M are well described for the p T > 25 
GeV by the sum of the LO contribution (J5]) and NLO contribution ([7]) after the mMRK 
subtraction. For the px > M data are well described by our prediction for all values of pr, 
and the NLO contribution is, in fact, negligible. The contribution of the box subprocess 


(17) is only about 15% of the cross section predicted at small pr, and decreases with px very 
fast, contributing significantly only for the px < 30 GeV. 

For the px < 25 GeV one can observe the deficit of the predicted cross section which 
reaches up to a factor of 5 at the px — 7 GeV. The region of small-p^ corresponds to the 
kinematics of CPM, where the radiation of soft gluons and virtual corrections are dominating. 
We expect, that computation of the NLO real-virtual interference correction in the PRA will 
significantly reduce this gap. One of the advantages of PRA is that at NLO this correction is 
finite and can be considered separately from the real NLO corrections, which are the subject 
of the present study. 

The good description of the data for the p T > M region supports the self-consistency 
of our approach, since the NLO correction in this region is almost canceled by the mMRK 
subtraction terms and the contribution of the real-virtual NLO correction is expected to be 
small here. 

For the reader’s convinience, in the figs. [8]-[T2j we also have plotted the corresponding 
NLO CPM predictions. Data for these plots correspond to the Diphox predictions HU, 
presented in the CDF ||3j and ATLAS (4] experimental papers. The contribution of the 
gg y 77 subprocess is also included into this predictions, via GAMMA2MC program m- 
Comparing the NLO CPM and NLO* PRA predictions in the figs. [8]- 12 one can conclude, 
that, NLO* approximation in PRA can not describe the data on the da/dM-distribution 
due to the absence of the loop correction, which contributes mostly in the back-to-back 
CPM-like kinematics. But for the configurations far away from the CPM kinematics, NLO* 
PRA describes data substantially better than NLO CPM, especially at the LffC. Moreover, 
in this region, the NLO* PRA prediction is dominated by the LO term, which demonstrates 
the better stability of PRA predictions for the kinematics far away from CPM one. Inclusion 
of full NLO corrections should also improve the agreement in the CPM region. 

In the left panel of the Fig. [9] one can observe the same qualitative features as in the 
left panel of the Fig. [8j despite the fact, that we have moved from Tevatron to the LHC 
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with it’s 3.6 times larger energy, and switched to pp collisions instead of pp ones. The NLO 
subprocess ([7]) is more important at the LHC than at the Tevatron, contributing significantly 
up to pt = 200 GeV. 

In the Fig. 10 and right panel of the Fig. [9] the A 0-spectra for the Tevatron and LHC 
are presented. In the both figures one can observe a good agreement of our predictions with 
data for A0 <1.5 which corresponds to the high deviation from the back-to-back kinematics 
for the photons. In this region the NLO correction is manifestly subleading, as it was for 
the pT-spectrum. The good description of the Tevatron data for the p T > M case is also 
there, as well, as the deficit of the predicted cross section for the back-to-back kinematics. 


As for the M-spectra of the Figs. 11 and 12 one certainly expects the deficit of the 
calculated cross section for the most values of M due to the deficit of the cross section for 
the CPM kinematics, observed earlier, since most of the total cross section is accumulated 
near the CPM configurations. However, in the region of M below the peak, the data are 
well described, demonstrating that PRA is suitable for the description of the effects of 
kinematical cuts. Once again, we observe the good description of the M-spectrum for the 
Pt > M subset of the Tevatron data. 


The contribution of the quark-box subprocess [17] to the M-spectra is found to be only 
about 8 % of the observed cross section in the peak, and 18% of the predicted cross section, 
both for the CDF-2012 and ATLAS-2013 kinematics. This result is 20—30% smaller than the 
usual CPM estimate [T2j, which is in accordance with the findings of the Ref. [52] where it 
was shown, that the space-like virtuality of the initial-state partons suppresses the 7 A —y yg 
contribution with respect to CPM expectation. 


VI. CONCLUSIONS 


In the present study, the pair hadroproduction of prompt photons is considered in the 
framework of PRA with tree-level NLO corrections (TtI | 8 l) and NNLO quark-box subprocess 


(17) taken into account. The procedure of localization in rapidity of the tree-level NLO 


corrections to avoid the double counting the real emissions between the hard-scattering part 


of the cross section and unPDF is proposed in the Sec. Ill As a consequence of this pro¬ 
cedure, the real NLO corrections are put under quantitative control, and their contribution 
was found to be numerically small at high px or in the kinematical region px > M. The 
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kinematical region pr > M is interesting for the further theoretical and experimental study, 
as an ideal testing site for the PRA, where the MRK between the 1SR and the hard sub¬ 
process is dominating. The contribution of the quark-box subprocess 0 was found to be 
about 8% of the observed cross section in the peak of the da/dM distribution, which is a 
bit smaller than the CPM estimate ra due to the space-like virtuality of the initial-state 
partons, similarly to the results of the Ref. [52]. 
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FIG. 2: 


The set of Feynman diagrams for the LO subprocess ([5]). 



FIG. 3: The set of the Feynman diagrams for the NLO subprocess Q. 
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FIG. 4: The set of the Feynman diagrams for the NLO subprocess Q. 




FIG. 5: Upper panel, the diagrammatic representation for the squared amplitude of the sub¬ 
process Q and the corresponding mMRK subtraction term. Lower panel - the same for the 
subprocess Q. 
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FIG. 6 : The Feynman diagrams for the mMRK subprodesses g(qi) + Q{q 2 ) ~> 9(93) + 9(94) (left 
panel) and 9(91) + Q{q 2 ) ->• 9(93) + 9(94) (right panel). 



FIG. 7 : The comparison of the px (left column) and Ay (right column) spectra fpr the NLO QMRK 
contributions (solid lines) with the corresponding mMRK subtraction terms (dashed lines), upper 
panels - subprocess 0 , lower panels - subprocess ([ 8 ]). 
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FIG. 8 : The py-spectra for the CDF -2012 dataset. The thick solid curve is the sum of the con¬ 
tributions ©> 0 with mMRK subtraction, and ©■ The thick dashed curve is the sum of the 
first two. The thin dash-dotted curve is the contribution of the subprocess ([ 5 ]) only. The thin 
dash-double-dotted curve is the corresponding Diphox (NLO CPM) prediction, taken from the 
Ref. [ 3 ], see the text for details. 

















































FIG. 9: The py-spectra (left panel) and A^-spectra (right panel) for the ATLAS-2013 dataset. 
The notation for the curves is the same as in the Fig. [8j except the thin dash-double-dotted curve, 
which is the Diphox (NLO CPM) prediction, taken from the Ref. @j, see the text for details. 
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FIG. 10: The A</>-spectra for the CDF-2012 dataset. The notation for the curves is the same as in 
the Fig. [8| 
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FIG. 11: Diphoton invariant mass spectra for the CDF-2012 dataset. The notation for the curves 
is the same as in the Fig. [8] 
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FIG. 12: Diphoton invariant mass spectrum for the ATLAS-2013 dataset. The notation for the 
curves is the same as in the Fig. [8] 
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